function [rmin,rmax] = genBounds(Q0,jshock,ivar,vma,N,Aineq,bineq,Aeq,...
    beq,LB,UB,optimOptions)

rmin = zeros(length(ivar),size(Q0,3));
rmax = rmin;

for kk = 1:size(Q0,3) % For each horizon
    
    QInit = Q0(:,:,kk);
    QInit = QInit(:); % Vectorise initial value of Q

    for ii = 1:length(ivar) % For each variable of interest

        % Minimise IRF subject to sign and equality restrictions.
        [~,rmin(ii,kk)] = fmincon(@(Q) genIRF(Q,jshock,...
            vma(ivar(ii),:),N,1),QInit,Aineq,bineq,Aeq,beq,LB,UB,...
            @(Q) Qcon(Q,N),optimOptions);
        % Maximise IRF subject to sign and equality restrictions.
        [~,rmax(ii,kk)] = fmincon(@(Q) genIRF(Q,jshock,...
            vma(ivar(ii),:),N,-1),QInit,Aineq,bineq,Aeq,beq,LB,UB,...
            @(Q) Qcon(Q,N),optimOptions);

    end

end